Run AlphaDiversity in scratchnotebooks
#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5 ✓ purrr 0.3.4
✓ tibble 3.1.6 ✓ dplyr 1.0.7
✓ tidyr 1.1.4 ✓ stringr 1.4.0
✓ readr 2.0.0 ✓ forcats 0.5.1
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
ksource(here::here("ActiveNotebooks", "AlphaDiversity.Rmd"))
processing file: /home/jacob/Projects/ChesapeaakeMainstemAnalysis/ActiveNotebooks/AlphaDiversity.Rmd
|
| | 0%
|
|...... | 3%
|
|............ | 7%
|
|.................. | 10%
|
|........................ | 13%
|
|.............................. | 17%
|
|................................... | 20%
|
|......................................... | 23%
|
|............................................... | 27%
|
|..................................................... | 30%
|
|........................................................... | 33%
|
|................................................................. | 37%
|
|....................................................................... | 40%
|
|............................................................................. | 43%
|
|................................................................................... | 47%
|
|........................................................................................ | 50%
|
|.............................................................................................. | 53%
|
|.................................................................................................... | 57%
|
|.......................................................................................................... | 60%
|
|................................................................................................................ | 63%
|
|...................................................................................................................... | 67%
|
|............................................................................................................................ | 70%
|
|.................................................................................................................................. | 73%
|
|........................................................................................................................................ | 77%
|
|.............................................................................................................................................. | 80%
|
|.................................................................................................................................................... | 83%
|
|......................................................................................................................................................... | 87%
|
|............................................................................................................................................................... | 90%
|
|..................................................................................................................................................................... | 93%
|
|........................................................................................................................................................................... | 97%
|
|.................................................................................................................................................................................| 100%
output file: /tmp/RtmpWB71Zc/file512e51983b76
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Attaching package: ‘flextable’
The following object is masked from ‘package:purrr’:
compose
Warning: Assuming taxa are rows
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3
library(cowplot)
library(flextable)
library(ftExtra)
This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.
Reshape back into an ASV matrix, but this time correcting for total abundance
nonSpikes %>% head
raDf <- nonSpikes %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")
raMat1 <- raMat %>% as.matrix()
countMat <- nonSpikes %>%
pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
column_to_rownames("ID") %>% as.matrix()
seqDep <- countMat %>% apply(1, sum)
min(seqDep)
[1] 852
sampleRichness <- rarefy(countMat, min(seqDep))
rarefy everything to the minimum depth (852)
countRare <- rrarefy(countMat, min(seqDep))
Gamma diversity
specpool(countRare)
Doesn’t finish
#specpool(countMat)
All richness estimates
richnessRare <- estimateR(countRare)
Shannon diversity
shan <- diversity(countRare)
shan
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2 3-1-S-180 3-1-S-20 3-1-S-5 3-2-B-0-2 3-2-B-1-2 3-2-B-180 3-2-B-20 3-2-B-5
4.266589 5.066675 4.748195 5.882121 5.213130 3.866813 5.519614 4.504019 4.997699 4.689925 5.265848 4.805358 4.315456 4.623881 4.559119 5.040415 5.448241
3-2-B-500 3-2-B-53 3-2-S-0-2 3-2-S-1-2 3-2-S-180 3-2-S-20 3-2-S-5 3-2-S-500 3-2-S-53 3-3-B-0-2 3-3-B-1-2 3-3-B-180 3-3-B-20 3-3-B-5 3-3-B-500 3-3-B-53 3-3-S-180
4.970669 4.754482 3.681555 4.892022 4.667467 5.000035 5.274235 4.862903 4.293725 4.391800 4.959762 3.460213 5.726654 5.328888 5.316755 5.443589 5.054444
3-3-S-20 3-3-S-500 3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180 4-3-B-20 4-3-B-5 4-3-B-500 4-3-B-53 4-3-O-1-2 4-3-O-180 4-3-O-5 4-3-O-500 4-3-O-53 4-3-S-0-2 4-3-S-180
4.906939 4.885977 4.377182 4.224096 4.916432 4.431618 4.408475 4.487108 4.300986 4.176164 4.968473 4.647133 5.112024 3.975860 4.573170 2.805419 4.498084
4-3-S-20 4-3-S-500 4-3-S-53 5-1-S-1-2 5-1-S-180 5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53 5-5-B-0-2 5-5-B-180 5-5-B-5 5-5-B-500 5-5-B-53 5-5-S-180 5-5-S-5 5-5-S-500
4.633447 4.763641 4.366296 4.403596 4.477610 4.171291 3.909744 4.041151 4.134044 4.665525 5.098082 5.398357 5.076910 5.036173 4.295120 4.994707 4.830118
5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2 C_5P1B_20 C_5P1B_500 C_5P1B_53
4.247241 4.407919 4.839992 4.888112 5.396460 4.679991 5.232058
Evenness
pielouJ <- shan/richnessRare["S.chao1",]
pielouJ
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2 3-1-S-180 3-1-S-20 3-1-S-5 3-2-B-0-2 3-2-B-1-2 3-2-B-180
0.010796521 0.008379501 0.012230611 0.003107806 0.006670746 0.042030574 0.008948671 0.011785878 0.006815941 0.008871841 0.007074017 0.007942428 0.013081311 0.011009240 0.012591089
3-2-B-20 3-2-B-5 3-2-B-500 3-2-B-53 3-2-S-0-2 3-2-S-1-2 3-2-S-180 3-2-S-20 3-2-S-5 3-2-S-500 3-2-S-53 3-3-B-0-2 3-3-B-1-2 3-3-B-180 3-3-B-20
0.006497564 0.006502422 0.007438591 0.008844683 0.019571778 0.006893855 0.008818684 0.008837454 0.009265714 0.009464584 0.014658199 0.014525960 0.006912331 0.062066595 0.005167280
3-3-B-5 3-3-B-500 3-3-B-53 3-3-S-180 3-3-S-20 3-3-S-500 3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180 4-3-B-20 4-3-B-5 4-3-B-500 4-3-B-53 4-3-O-1-2
0.006820821 0.007008356 0.006606406 0.010298444 0.010244131 0.010936137 0.011279040 0.012088416 0.009557943 0.011023255 0.009479392 0.008788068 0.010645362 0.009541883 0.008057076
4-3-O-180 4-3-O-5 4-3-O-500 4-3-O-53 4-3-S-0-2 4-3-S-180 4-3-S-20 4-3-S-500 4-3-S-53 5-1-S-1-2 5-1-S-180 5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53
0.010062029 0.008070617 0.015028669 0.009695583 0.124685269 0.010941017 0.008496069 0.010528326 0.013690750 0.010095361 0.012491681 0.013629647 0.017070925 0.014583301 0.015952761
5-5-B-0-2 5-5-B-180 5-5-B-5 5-5-B-500 5-5-B-53 5-5-S-180 5-5-S-5 5-5-S-500 5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2 C_5P1B_20 C_5P1B_500 C_5P1B_53
0.008436353 0.007538475 0.006170381 0.011690997 0.009564602 0.014591206 0.009765933 0.009603814 0.014900732 0.006386277 0.007006523 0.006184548 0.003893466 0.008552268 0.005238606
diversityData <- sampleData %>%
left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
arrange(Size_Class)
Parameters for all plots
plotSpecs <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
Observed species counts, on rarefied data
plotObs <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs
Estemated Chao1 Richness
plotChao1 <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) +
scale_y_log10() +
ylab("Richness (Chao1)")
plotChao1
Shannon diversity
plotShan <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
ylab("Diversity (Shannon H)") +
lims(y = c(2.5, 6))
plotShan
Evenness
plotPielou <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou
All plots together
plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha
ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)
Rarefied
obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)
Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-222.130 -40.254 0.968 47.519 200.464
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 294642.461 134921.667 2.184 0.032379 *
log(Size_Class) 33.747 8.174 4.128 0.000101 ***
I(log(Size_Class)^2) -6.605 1.521 -4.343 4.72e-05 ***
lat -15345.129 7030.594 -2.183 0.032470 *
I(lat^2) 199.852 91.525 2.184 0.032397 *
depth 5.038 3.238 1.556 0.124293
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 76.85 on 69 degrees of freedom
Multiple R-squared: 0.2783, Adjusted R-squared: 0.226
F-statistic: 5.322 on 5 and 69 DF, p-value: 0.0003409
chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-516.16 -139.22 -14.04 140.31 1146.14
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 952088.842 448528.950 2.123 0.037370 *
log(Size_Class) 85.540 27.175 3.148 0.002430 **
I(log(Size_Class)^2) -18.549 5.056 -3.669 0.000476 ***
lat -49615.130 23372.263 -2.123 0.037359 *
I(lat^2) 646.322 304.262 2.124 0.037237 *
depth 19.787 10.763 1.838 0.070303 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 255.5 on 69 degrees of freedom
Multiple R-squared: 0.2193, Adjusted R-squared: 0.1627
F-statistic: 3.877 on 5 and 69 DF, p-value: 0.003731
chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2),
data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-451.11 -148.26 -32.91 126.61 1234.57
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 569.961 46.362 12.294 < 2e-16 ***
log(Size_Class) 86.118 27.533 3.128 0.002542 **
I(log(Size_Class)^2) -18.923 5.122 -3.695 0.000426 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 259.2 on 72 degrees of freedom
Multiple R-squared: 0.1617, Adjusted R-squared: 0.1384
F-statistic: 6.943 on 2 and 72 DF, p-value: 0.001748
shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(shanMod)
Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-1.32009 -0.22746 0.04123 0.30999 0.76552
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.887e+03 7.925e+02 2.382 0.0200 *
log(Size_Class) 2.051e-01 4.801e-02 4.273 6.06e-05 ***
I(log(Size_Class)^2) -3.764e-02 8.933e-03 -4.213 7.48e-05 ***
lat -9.807e+01 4.129e+01 -2.375 0.0203 *
I(lat^2) 1.277e+00 5.376e-01 2.375 0.0204 *
depth 2.466e-02 1.902e-02 1.297 0.1991
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4514 on 69 degrees of freedom
Multiple R-squared: 0.3087, Adjusted R-squared: 0.2586
F-statistic: 6.161 on 5 and 69 DF, p-value: 8.916e-05
pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)
Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-0.016223 -0.005097 -0.002264 0.000799 0.099586
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.365e+01 2.621e+01 -0.902 0.3700
log(Size_Class) -4.048e-03 1.588e-03 -2.549 0.0130 *
I(log(Size_Class)^2) 7.103e-04 2.955e-04 2.404 0.0189 *
lat 1.232e+00 1.366e+00 0.902 0.3704
I(lat^2) -1.601e-02 1.778e-02 -0.901 0.3709
depth -3.282e-04 6.290e-04 -0.522 0.6034
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01493 on 69 degrees of freedom
Multiple R-squared: 0.1032, Adjusted R-squared: 0.03826
F-statistic: 1.589 on 5 and 69 DF, p-value: 0.1748
uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2
predict(obsMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
225.2948 225.2948 191.7909 191.7909 205.5804 152.9884 152.9884 184.9565 198.0723 302.6501 302.6501 269.1462 269.1462 282.9357 230.3437 230.3437 262.3118 262.3118 333.9210 333.9210
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
300.4171 300.4171 314.2066 261.6146 261.6146 293.5827 306.6985 306.6985 338.5363 338.5363 305.0324 305.0324 318.8220 318.8220 266.2299 266.2299 298.1980 298.1980 326.5841 293.0801
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
293.0801 306.8697 306.8697 254.2777 254.2777 254.2777 286.2457 286.2457 299.3615 299.3615 293.8447 293.8447 260.3408 260.3408 274.1303 274.1303 221.5383 221.5383 221.5383 253.5064
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
253.5064 266.6221 266.6221 251.3431 217.8392 217.8392 231.6287 231.6287 179.0367 179.0367 179.0367 211.0048 211.0048 224.1206 224.1206
$se.fit
[1] 26.88498 26.88498 26.15198 26.15198 26.01242 27.67418 27.67418 29.22809 33.90698 18.85945 18.85945 18.24545 18.24545 17.13460 19.96313 19.96313 21.42036 21.42036 19.17995 19.17995
[21] 18.71991 18.71991 17.16981 20.00965 20.00965 21.12357 28.06863 28.06863 19.17779 19.17779 18.68694 18.68694 16.89659 16.89659 19.53551 19.53551 20.49627 20.49627 18.49063 17.85404
[41] 17.85404 15.93670 15.93670 18.37221 18.37221 18.37221 19.35941 19.35941 26.45249 26.45249 19.23795 19.23795 18.35438 18.35438 16.61635 16.61635 18.35588 18.35588 18.35588 19.42396
[61] 19.42396 26.01457 26.01457 24.36577 23.41633 23.41633 22.25921 22.25921 23.05203 23.05203 23.05203 24.04762 24.04762 29.11453 29.11453
$df
[1] 69
$residual.scale
[1] 76.8497
diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit
plotSpecs2 <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
#geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
plotObs_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted ASVs")
plotObs_pred
predict(chao1Mod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
471.0415 471.0415 357.8857 357.8857 442.2027 277.6790 277.6790 403.5903 380.9771 671.7413 671.7413 558.5855 558.5855 642.9025 478.3788 478.3788 604.2901 604.2901 746.3854 746.3854
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
633.2296 633.2296 717.5466 553.0229 553.0229 678.9342 656.3210 656.3210 746.5474 746.5474 633.3915 633.3915 717.7085 717.7085 553.1848 553.1848 679.0961 679.0961 703.9826 590.8268
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
590.8268 675.1438 675.1438 510.6201 510.6201 510.6201 636.5314 636.5314 613.9182 613.9182 600.7490 600.7490 487.5932 487.5932 571.9102 571.9102 407.3865 407.3865 407.3865 533.2978
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
533.2978 510.6846 510.6846 471.9554 358.7995 358.7995 443.1165 443.1165 278.5929 278.5929 278.5929 404.5041 404.5041 381.8910 381.8910
$se.fit
[1] 89.37550 89.37550 86.93875 86.93875 86.47479 91.99908 91.99908 97.16485 112.71919 62.69571 62.69571 60.65454 60.65454 56.96169 66.36473 66.36473 71.20912 71.20912
[19] 63.76118 63.76118 62.23181 62.23181 57.07873 66.51939 66.51939 70.22248 93.31038 93.31038 63.75398 63.75398 62.12222 62.12222 56.17044 56.17044 64.94318 64.94318
[37] 68.13710 68.13710 61.46961 59.35336 59.35336 52.97941 52.97941 61.07595 61.07595 61.07595 64.35777 64.35777 87.93774 87.93774 63.95397 63.95397 61.01666 61.01666
[55] 55.23882 55.23882 61.02166 61.02166 61.02166 64.57236 64.57236 86.48194 86.48194 81.00074 77.84443 77.84443 73.99777 73.99777 76.63338 76.63338 76.63338 79.94309
[73] 79.94309 96.78733 96.78733
$df
[1] 69
$residual.scale
[1] 255.4765
diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit
plotChao1_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred
predict(shanMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
4.495853 4.495853 4.294730 4.294730 4.302302 3.967491 3.967491 4.130228 4.372581 4.959672 4.959672 4.758549 4.758549 4.766122 4.431310 4.431310 4.594047 4.594047 5.156199 5.156199
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
4.955075 4.955075 4.962648 4.627837 4.627837 4.790574 5.032926 5.032926 5.200306 5.200306 4.999182 4.999182 5.006755 5.006755 4.671944 4.671944 4.834681 4.834681 5.144714 4.943591
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
4.943591 4.951164 4.951164 4.616353 4.616353 4.616353 4.779090 4.779090 5.021442 5.021442 4.973856 4.973856 4.772732 4.772732 4.780305 4.780305 4.445494 4.445494 4.445494 4.608231
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
4.608231 4.850584 4.850584 4.744786 4.543663 4.543663 4.551236 4.551236 4.216425 4.216425 4.216425 4.379161 4.379161 4.621514 4.621514
$se.fit
[1] 0.15791178 0.15791178 0.15360644 0.15360644 0.15278670 0.16254721 0.16254721 0.17167427 0.19915623 0.11077297 0.11077297 0.10716658 0.10716658 0.10064191 0.11725554 0.11725554
[17] 0.12581478 0.12581478 0.11265549 0.11265549 0.10995335 0.10995335 0.10084870 0.11752879 0.11752879 0.12407155 0.16486406 0.16486406 0.11264276 0.11264276 0.10975972 0.10975972
[33] 0.09924390 0.09924390 0.11474388 0.11474388 0.12038702 0.12038702 0.10860666 0.10486760 0.10486760 0.09360588 0.09360588 0.10791114 0.10791114 0.10791114 0.11370956 0.11370956
[49] 0.15537149 0.15537149 0.11299611 0.11299611 0.10780639 0.10780639 0.09759788 0.09759788 0.10781521 0.10781521 0.10781521 0.11408871 0.11408871 0.15279933 0.15279933 0.14311495
[65] 0.13753827 0.13753827 0.13074185 0.13074185 0.13539854 0.13539854 0.13539854 0.14124626 0.14124626 0.17100726 0.17100726
$df
[1] 69
$residual.scale
[1] 0.4513849
diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit
plotShannonH_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred
predict(pielouMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
0.019650593 0.019650593 0.022081017 0.022081017 0.021538847 0.025099422 0.025099422 0.022609276 0.019055594 0.010581192 0.010581192 0.013011616 0.013011616 0.012469446 0.016030021
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
0.016030021 0.013539875 0.013539875 0.006620266 0.006620266 0.009050690 0.009050690 0.008508520 0.012069095 0.012069095 0.009578950 0.006025267 0.006025267 0.005542827 0.005542827
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
0.007973250 0.007973250 0.007431081 0.007431081 0.010991656 0.010991656 0.008501510 0.008501510 0.006419563 0.008849987 0.008849987 0.008307817 0.008307817 0.011868392 0.011868392
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
0.011868392 0.009378246 0.009378246 0.005824564 0.005824564 0.009427605 0.009427605 0.011858029 0.011858029 0.011315859 0.011315859 0.014876434 0.014876434 0.014876434 0.012386288
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
0.012386288 0.008832606 0.008832606 0.013569678 0.016000102 0.016000102 0.015457932 0.015457932 0.019018507 0.019018507 0.019018507 0.016528362 0.016528362 0.012974679 0.012974679
$se.fit
[1] 0.005223138 0.005223138 0.005080733 0.005080733 0.005053619 0.005376461 0.005376461 0.005678351 0.006587352 0.003663961 0.003663961 0.003544674 0.003544674 0.003328862 0.003878380
[16] 0.003878380 0.004161488 0.004161488 0.003726227 0.003726227 0.003636850 0.003636850 0.003335702 0.003887418 0.003887418 0.004103828 0.005453094 0.005453094 0.003725806 0.003725806
[31] 0.003630446 0.003630446 0.003282621 0.003282621 0.003795304 0.003795304 0.003981958 0.003981958 0.003592307 0.003468633 0.003468633 0.003096136 0.003096136 0.003569302 0.003569302
[46] 0.003569302 0.003761092 0.003761092 0.005139115 0.005139115 0.003737494 0.003737494 0.003565837 0.003565837 0.003228177 0.003228177 0.003566129 0.003566129 0.003566129 0.003773633
[61] 0.003773633 0.005054037 0.005054037 0.004733714 0.004549258 0.004549258 0.004324457 0.004324457 0.004478483 0.004478483 0.004478483 0.004671904 0.004671904 0.005656288 0.005656288
$df
[1] 69
$residual.scale
[1] 0.01493014
diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit
plot_pielouJ_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred
plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).
plotPredictions
ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)
plot_grid(plotObs, plotChao1, plotShan, plotPielou,
plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
nrow = 2, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).
alphaSummary <- tibble(
metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(obsMod, chao1Mod, shanMod, pielouMod)
)
alphaSummary <- alphaSummary %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
alphaSummary <- alphaSummary %>%
select(-model) %>%
unnest(df)
alphaSummary <- alphaSummary %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
Metric | Term | Estimate | Standard | T Value | p |
Observed ASVs | Intercept | 2.9 x 105 | 1.3 x 105 | 2.18 | 0.032 |
log(Size Class) | 3.4 x 101 | 8.2 x 100 | 4.13 | < 0.001 | |
log(Size Class)2 | -6.6 x 100 | 1.5 x 100 | -4.34 | < 0.001 | |
Latitude | -1.5 x 104 | 7.0 x 103 | -2.18 | 0.032 | |
Latitude2 | 2.0 x 102 | 9.2 x 101 | 2.18 | 0.032 | |
Depth | 5.0 x 100 | 3.2 x 100 | 1.56 | 0.124 | |
Richness (Chao1) | Intercept | 9.5 x 105 | 4.5 x 105 | 2.12 | 0.037 |
log(Size Class) | 8.6 x 101 | 2.7 x 101 | 3.15 | 0.002 | |
log(Size Class)2 | -1.9 x 101 | 5.1 x 100 | -3.67 | < 0.001 | |
Latitude | -5.0 x 104 | 2.3 x 104 | -2.12 | 0.037 | |
Latitude2 | 6.5 x 102 | 3.0 x 102 | 2.12 | 0.037 | |
Depth | 2.0 x 101 | 1.1 x 101 | 1.84 | 0.070 | |
Diversity (Shannon H) | Intercept | 1.9 x 103 | 7.9 x 102 | 2.38 | 0.020 |
log(Size Class) | 2.1 x 10-1 | 4.8 x 10-2 | 4.27 | < 0.001 | |
log(Size Class)2 | -3.8 x 10-2 | 8.9 x 10-3 | -4.21 | < 0.001 | |
Latitude | -9.8 x 101 | 4.1 x 101 | -2.37 | 0.020 | |
Latitude2 | 1.3 x 100 | 5.4 x 10-1 | 2.37 | 0.020 | |
Depth | 2.5 x 10-2 | 1.9 x 10-2 | 1.30 | 0.199 | |
Evenness (Pielou J) | Intercept | -2.4 x 101 | 2.6 x 101 | -0.90 | 0.370 |
log(Size Class) | -4.0 x 10-3 | 1.6 x 10-3 | -2.55 | 0.013 | |
log(Size Class)2 | 7.1 x 10-4 | 3.0 x 10-4 | 2.40 | 0.019 | |
Latitude | 1.2 x 100 | 1.4 x 100 | 0.90 | 0.370 | |
Latitude2 | -1.6 x 10-2 | 1.8 x 10-2 | -0.90 | 0.371 | |
Depth | -3.3 x 10-4 | 6.3 x 10-4 | -0.52 | 0.603 |
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))
Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.
diversityDataWB <- full_join(diversityData,
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
mutate(pielouJ2 = shannonH/break_estimate) %>%
identity()
diversityDataWB
pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)
Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-0.013948 -0.005120 -0.002519 0.000893 0.105997
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.940e+01 2.756e+01 -0.704 0.4840
log(Size_Class) -3.284e-03 1.670e-03 -1.966 0.0533 .
I(log(Size_Class)^2) 5.739e-04 3.107e-04 1.847 0.0690 .
lat 1.009e+00 1.436e+00 0.702 0.4849
I(lat^2) -1.310e-02 1.870e-02 -0.701 0.4858
depth -2.331e-04 6.614e-04 -0.352 0.7256
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.0157 on 69 degrees of freedom
Multiple R-squared: 0.06678, Adjusted R-squared: -0.0008485
F-statistic: 0.9875 on 5 and 69 DF, p-value: 0.4318
Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.
plotBreak <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Richness (Breakaway)")
plotBreak
plotPielou2 <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Evenness (PielouJ)")
plotPielou2
predict(pielouMod2, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13
0.0116712168 0.0116712168 0.0135414829 0.0135414829 0.0133652561 0.0159774714 0.0159774714 0.0139770206 0.0099002360 0.0043201718 0.0043201718 0.0061904378 0.0061904378
14 15 16 17 18 19 20 21 22 23 24 25 26
0.0060142111 0.0086264264 0.0086264264 0.0066259755 0.0066259755 0.0011013110 0.0011013110 0.0029715770 0.0029715770 0.0027953503 0.0054075656 0.0054075656 0.0034071147
27 28 29 30 31 32 33 34 35 36 37 38 39
-0.0006696698 -0.0006696698 0.0002127168 0.0002127168 0.0020829828 0.0020829828 0.0019067561 0.0019067561 0.0045189714 0.0045189714 0.0025185205 0.0025185205 0.0009083804
40 41 42 43 44 45 46 47 48 49 50 51 52
0.0027786464 0.0027786464 0.0026024197 0.0026024197 0.0052146350 0.0052146350 0.0052146350 0.0032141841 0.0032141841 -0.0008626004 -0.0008626004 0.0033228151 0.0033228151
53 54 55 56 57 58 59 60 61 62 63 64 65
0.0051930811 0.0051930811 0.0050168543 0.0050168543 0.0076290696 0.0076290696 0.0076290696 0.0056286188 0.0056286188 0.0015518342 0.0015518342 0.0066561189 0.0085263849
66 67 68 69 70 71 72 73 74 75
0.0085263849 0.0083501582 0.0083501582 0.0109623735 0.0109623735 0.0109623735 0.0089619227 0.0089619227 0.0048851381 0.0048851381
$se.fit
[1] 0.005492160 0.005492160 0.005342421 0.005342421 0.005313910 0.005653380 0.005653380 0.005970819 0.006926639 0.003852676 0.003852676 0.003727246 0.003727246 0.003500318 0.004078139
[16] 0.004078139 0.004375829 0.004375829 0.003918150 0.003918150 0.003824170 0.003824170 0.003507510 0.004087643 0.004087643 0.004315200 0.005733960 0.005733960 0.003917707 0.003917707
[31] 0.003817435 0.003817435 0.003451696 0.003451696 0.003990784 0.003990784 0.004187052 0.004187052 0.003777332 0.003647288 0.003647288 0.003255606 0.003255606 0.003753142 0.003753142
[46] 0.003753142 0.003954810 0.003954810 0.005403809 0.005403809 0.003929997 0.003929997 0.003749498 0.003749498 0.003394447 0.003394447 0.003749805 0.003749805 0.003749805 0.003967997
[61] 0.003967997 0.005314350 0.005314350 0.004977528 0.004783571 0.004783571 0.004547192 0.004547192 0.004709152 0.004709152 0.004709152 0.004912535 0.004912535 0.005947620 0.005947620
$df
[1] 69
$residual.scale
[1] 0.01569913
diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit
plot_pielouJ2_pred <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred
plotBreakaway <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) +
scale_y_log10() +
ylab("Richness (Breakaway)")
plotBreakaway
#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame
Why are these not smooth curves?!! What if I redo the model, this time with the same data frame
breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat + I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()
Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-2974.5 -1191.2 -151.6 599.9 6768.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7124615.61 3339862.88 2.133 0.0365 *
log(Size_Class) 244.45 202.35 1.208 0.2312
I(log(Size_Class)^2) -75.16 37.65 -1.996 0.0498 *
lat -370568.38 174035.93 -2.129 0.0368 *
I(lat^2) 4817.28 2265.61 2.126 0.0371 *
depth 151.10 80.15 1.885 0.0636 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared: 0.1414, Adjusted R-squared: 0.0792
F-statistic: 2.273 on 5 and 69 DF, p-value: 0.0567
Note the non statistical significance overall
#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit
plot_break_pred <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
# filter(Station == 4.3) %>%
ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred
plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB
ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlpha, width = 11, height = 4)
Summary table I want both breakaway metrics here
bettaTable <- myBet$table %>%
as.data.frame() %>%
rename(estimate = Estimates,
`std.error` = `Standard Errors`,
`p.value`=`p-values`
) %>%
mutate(`statistic` = NA) %>%
rownames_to_column(var = "term") %>%
select(term, estimate, std.error, statistic, p.value) %>%
as_tibble()
bettaTable
alphaSummary2 <- tibble(
metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(breakLm, shanMod, pielouMod2)
)
alphaSummary2 <- alphaSummary2 %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
## Add in willis variables
breakawaySummary <- tibble(
metric = "Richness (Breakaway -- Betta)",
model = NULL,
df = list(bettaTable)
)
alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)
alphaSummary2 <- alphaSummary2 %>%
select(-model) %>%
unnest(df)
alphaSummary2 <- alphaSummary2 %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary2
alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
alphaTable2
Metric | Term | Estimate | Standard | T Value | p |
Richness (Breakaway – Betta) | Intercept | 7.1 x 106 | 2.4 x 102 | NA | < 0.001 |
log(Size Class) | 1.2 x 102 | 6.1 x 101 | NA | 0.058 | |
log(Size Class)2 | -5.0 x 101 | 1.2 x 101 | NA | < 0.001 | |
Latitude | -3.7 x 105 | 6.1 x 100 | NA | < 0.001 | |
Latitude2 | 4.8 x 103 | 1.6 x 10-1 | NA | < 0.001 | |
Depth | 1.5 x 102 | 1.0 x 101 | NA | < 0.001 | |
Richness (Breakaway – LM) | Intercept | 7.1 x 106 | 3.3 x 106 | 2.13 | 0.036 |
log(Size Class) | 2.4 x 102 | 2.0 x 102 | 1.21 | 0.231 | |
log(Size Class)2 | -7.5 x 101 | 3.8 x 101 | -2.00 | 0.050 | |
Latitude | -3.7 x 105 | 1.7 x 105 | -2.13 | 0.037 | |
Latitude2 | 4.8 x 103 | 2.3 x 103 | 2.13 | 0.037 | |
Depth | 1.5 x 102 | 8.0 x 101 | 1.89 | 0.064 | |
Diversity (Shannon H) | Intercept | 1.9 x 103 | 7.9 x 102 | 2.38 | 0.020 |
log(Size Class) | 2.1 x 10-1 | 4.8 x 10-2 | 4.27 | < 0.001 | |
log(Size Class)2 | -3.8 x 10-2 | 8.9 x 10-3 | -4.21 | < 0.001 | |
Latitude | -9.8 x 101 | 4.1 x 101 | -2.37 | 0.020 | |
Latitude2 | 1.3 x 100 | 5.4 x 10-1 | 2.37 | 0.020 | |
Depth | 2.5 x 10-2 | 1.9 x 10-2 | 1.30 | 0.199 | |
Evenness (Pielou J) | Intercept | -1.9 x 101 | 2.8 x 101 | -0.70 | 0.484 |
log(Size Class) | -3.3 x 10-3 | 1.7 x 10-3 | -1.97 | 0.053 | |
log(Size Class)2 | 5.7 x 10-4 | 3.1 x 10-4 | 1.85 | 0.069 | |
Latitude | 1.0 x 100 | 1.4 x 100 | 0.70 | 0.485 | |
Latitude2 | -1.3 x 10-2 | 1.9 x 10-2 | -0.70 | 0.486 | |
Depth | -2.3 x 10-4 | 6.6 x 10-4 | -0.35 | 0.726 |
alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))
myBet$table
plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)